home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / lu_complex.pro < prev    next >
Text File  |  1997-07-08  |  7KB  |  176 lines
  1. ;$Id: lu_complex.pro,v 1.9 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       LU_COMPLEX
  8. ;
  9. ; PURPOSE:
  10. ;       This function solves an N by N complex linear system using
  11. ;       LU decomposition. The result is an N-element complex vector.
  12. ;       Alternatively, this function computes the generalized inverse 
  13. ;       of an N by N complex array using LU decomposition. The result 
  14. ;       is an N by N complex array.
  15. ;
  16. ; CATEGORY:
  17. ;       Complex Linear Algebra.
  18. ;
  19. ; CALLING SEQUENCE:
  20. ;       Result = LU_COMPLEX(A, B)
  21. ;
  22. ; INPUTS:
  23. ;       A:    An N by N complex array.
  24. ;
  25. ;       B:    An N-element right-side vector (real or complex).
  26. ;     
  27. ; KEYWORD PARAMETERS:
  28. ;       DOUBLE: If set to a non-zero value, computations are done in
  29. ;               double precision arithmetic.
  30. ;
  31. ;      INVERSE: If set to a non-zero value, the generalized inverse of A
  32. ;               is computed. In this case the input parameter B is ignored.
  33. ;
  34. ;       SPARSE: If set to a non-zero value, the input array is converted
  35. ;               to row-indexed sparse storage format. Computations are
  36. ;               done using the iterative biconjugate gradient method.
  37. ;               This keyword is effective only when solving complex linear
  38. ;               systems. This keyword has no effect when calculating the
  39. ;               generalized inverse.
  40. ;
  41. ; EXAMPLE:
  42. ;       1) Define a complex array (A) and right-side vector (B).
  43. ;            A = [[complex(1, 0), complex(2,-2), complex(-3,1)], $
  44. ;                 [complex(1,-2), complex(2, 2), complex(1, 0)], $
  45. ;                 [complex(1, 1), complex(0, 1), complex(1, 5)]]
  46. ;            B =  [complex(1, 1), complex(3,-2), complex(1,-2)]
  47. ;
  48. ;          Solve the complex linear system (Az = B) for z.
  49. ;            z = LU_COMPLEX(a, b)
  50. ;
  51. ;        2) Compute the generalized inverse of A.
  52. ;            inv = LU_COMPLEX(a, b, /inverse)
  53. ;
  54. ; PROCEDURE:
  55. ;       LU_COMPLEX solves the complex linear system Az = b using
  56. ;       LU decomposition. If the SPARSE keyword is set, the coefficient
  57. ;       array is converted to row-indexed sparse storage format and the 
  58. ;       system is solved using the iterative biconjugate gradient method.
  59. ;       LU_COMPLEX computes the generalized inverse of the complex
  60. ;       array A using LU decomposition if B is supplied as an arbitrary
  61. ;       scalar value or if the INVERSE keyword is set.
  62. ;
  63. ; REFERENCE:
  64. ;       Numerical Recipes, The Art of Scientific Computing (Second Edition)
  65. ;       Cambridge University Press
  66. ;       ISBN 0-521-43108-5
  67. ;
  68. ; MODIFICATION HISTORY:
  69. ;       Written by:  GGS, RSI, October 1993
  70. ;       Modified:    GGS, RSI, February 1994
  71. ;                    Transposing the array prior to calling LU_COMPLEX
  72. ;                    is no longer necessary. LU_COMPLEX is now able to
  73. ;                    compute the generalized inverse of an N by N complex
  74. ;                    array using LU decomposition.
  75. ;       Modified:    GGS, RSI, June 1994
  76. ;                    Included support for sparse complex arrays using the
  77. ;                    Numerical Recipes functions NR_SPRSIN and NR_LINBCG.
  78. ;       Modified:    GGS, RSI, Decemberber 1994
  79. ;                    Added support for double-precision complex inputs.
  80. ;                    Reduced internal memory allocation requirements.
  81. ;                    Added INVERSE keyword. New documentation header.
  82. ;       Modified:    GGS, RSI, April 1996
  83. ;                    Modified keyword checking and use of double precision. 
  84. ;-
  85.  
  86. FUNCTION LU_Complex, A, B, Double = Double, Inverse = Inverse, Sparse = Sparse
  87.  
  88.   ON_ERROR, 2  ;Return to caller if error occurs.
  89.  
  90.   if N_PARAMS() ne 2 then $
  91.     MESSAGE, "Incorrect number of input arguments."
  92.  
  93.   TypeA = SIZE(A)
  94.   TypeB = SIZE(B)
  95.  
  96.   if TypeA[TypeA[0]+1] ne 6 and TypeA[TypeA[0]+1] ne 9 then $
  97.     MESSAGE, "Input array must be complex or double-complex."
  98.  
  99.   if TypeA[1] ne TypeA[2] then $
  100.     MESSAGE, "Input array must be square."
  101.  
  102.   if TypeA[2] ne TypeB[TypeB[0]+2] and TypeB[TypeB[0]+2] ne 1 then $
  103.     MESSAGE, "Input array and right-side vector are of incompatible size."
  104.  
  105.   ;If the DOUBLE keyword is not set then the internal precision and
  106.   ;result are determined by the type of input.
  107.   if N_ELEMENTS(Double) eq 0 then $
  108.     Double = (TypeA[TypeA[0]+1] eq 9 or TypeB[TypeB[0]+1] eq 5 or $
  109.                                         TypeB[TypeB[0]+1] eq 9)
  110.  
  111.   ;Double-precision complex.
  112.   if Double ne 0 then begin
  113.     Comp = [[DOUBLE(A), -IMAGINARY(A)], $
  114.             [IMAGINARY(A), DOUBLE(A)]]
  115.     ;Generalized inverse of A (does not depend upon SPARSE keyword).
  116.     if TypeB[TypeB[0]+2] eq 1 or KEYWORD_SET(Inverse) then begin
  117.       Vec = DBLARR(2L*TypeA[1])
  118.       Inv = DCOMPLEXARR(TypeA[1], TypeA[1])
  119.       ;Compute the LU decomp only once and iterate on it!
  120.       LUDC, Comp, Index, Double = Double
  121.       for k = 0, TypeA[1]-1 do begin
  122.         Vec[k] = 1.0d
  123.         Sol = LUSOL(Comp, Index, Vec, Double = Double)
  124.         Vec[k] = 0.0d
  125.         Inv[k, *] = DCOMPLEX(Sol[0:TypeA[1]-1], Sol[TypeA[1]:*])
  126.       endfor
  127.       RETURN, Inv
  128.     endif else begin ;Solve Az = b
  129.       ;Rhs complex?
  130.       if TypeB[TypeB[0]+1] ne 6 and TypeB[TypeB[0]+1] ne 9 then $
  131.         Vec = [B, DBLARR(TypeB[TypeB[0]+2])] $ ;No.
  132.       else Vec = [DOUBLE(B), IMAGINARY(B)] ;Yes.
  133.       if keyword_set(SPARSE) eq 0 then begin ;Dense coefficient array.
  134.         LUDC, Comp, Index, Double = Double
  135.         Sol = LUSOL(Comp, Index, Vec, Double = Double)
  136.       endif else begin ;Sparse coefficient array.
  137.         Sol = LINBCG(SPRSIN(Comp, Double = Double), Vec, $
  138.                      REPLICATE(1.0d, 2L*TypeB[TypeB[0]+2]), Double = Double)
  139.       endelse
  140.       RETURN, DCOMPLEX(Sol[0:TypeA[1]-1], Sol[TypeA[1]:*])
  141.     endelse
  142.   ;Single-precision complex.
  143.   endif else begin
  144.     Comp = [[FLOAT(a), -IMAGINARY(a)], $
  145.             [IMAGINARY(a), FLOAT(a)]]
  146.     ;Generalized Inverse of A (does not depend upon SPARSE keyword).
  147.     if TypeB[TypeB[0]+2] eq 1 or KEYWORD_SET(Inverse) then begin
  148.       Vec = fltarr(2L*TypeA[1])
  149.       Inv = complexarr(TypeA[1], TypeA[1])
  150.       ;Compute the LU decomp only once and iterate on it!
  151.       LUDC, Comp, Index, Double = Double
  152.       for k = 0, TypeA[1]-1 do begin
  153.         Vec[k] = 1.0
  154.         Sol = LUSOL(Comp, Index, Vec, Double = Double)
  155.         Vec[k] = 0.0
  156.         Inv[k, *] = COMPLEX(Sol[0:TypeA[1]-1], Sol[TypeA[1]:*])
  157.       endfor
  158.       RETURN, Inv
  159.     endif else begin ;Solve Az = b
  160.       ;Rhs complex?
  161.       if TypeB[TypeB[0]+1] ne 6 and TypeB[TypeB[0]+1] ne 9 then $
  162.         Vec = [B, FLTARR(TypeB[TypeB[0]+2])] $  ;No.
  163.         else Vec = [FLOAT(B), IMAGINARY(B)] ;Yes.
  164.       if keyword_set(SPARSE) eq 0 then begin ;Dense coefficient array.
  165.         LUDC, Comp, Index, Double = Double
  166.         Sol = LUSOL(Comp, Index, Vec, Double = Double)
  167.       endif else begin ;Sparse coefficient array.
  168.         Sol = LINBCG(SPRSIN(Comp, Double = Double), Vec, $
  169.                      REPLICATE(1.0, 2L*TypeB[TypeB[0]+2]), Double = Double)
  170.       endelse
  171.       RETURN, COMPLEX(Sol[0:TypeA[1]-1], Sol[TypeA[1]:*])
  172.     endelse
  173.   endelse
  174.  
  175. END
  176.